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Abstract. The quantum phase transition point between the insulator and the superfluid phase at unit 
filling factor of the infinite one-dimensional Bose-Hubbard model is numerically computed with a 
high accuracy. The method uses the infinite system version of the time evolving block decimation 
algorithm, here tested in a challenging case. We provide also the accurate estimate of the phase 
transition point at double occupancy. 
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INTRODUCTION 

Ultra-cold atoms in optical lattices make possible detailed studies of basic condensed 
matter systems [Ul2l[3l|4||. While several applications have been envisaged ranging from 
high T c superconductivity [5], disordered systems [6], spin models [7] or novel quantum 
magnets [8) - for a review, see - the most basic model realized in such a setting is 
the Bose-Hubbard model (BHM) [101, which still defies an exact analytical description 
while exhibiting a nontrivial quantum phase transition [fTTTl between a superfluid state 
and an insulator. Experiments, both three-dimensional and one-dimensional (ID) 
0, performed in additional confining atomic traps, led to realization of an inhomoge- 
neous finite model. Still, from the theoretical perspective, the infinite system is the most 
challenging case. As in experiments, the typical numerical studies deal with finite size 
systems, with a possible extrapolation to infinite systems at the end. Such an extrapo- 
lation is always delicate, especially in the vicinity of the phase transition. Even simple 
quantities like the precise location of the phase transition point are still controversial, 
see below. 

In the absence of an exact analytical result, one must rely on approximate results, such 
as the Bethe ansatz [f!2ll or numerical studies, which are roughly of three different types: 
exact diagonalization for finite systems - limited to very small systems and thus unable 
to extrapolate to the infinite size -, Quantum Monte-Carlo (QMC) methods and Density 
Matrix Renormalization Group (DMRG) methods lfT3i 

The Time Evolving Block Decimation (TEBD) algorithm [fT4l [T51l is a novel variant 
of the DMRG method which can provide results for the temporal dynamics, but also 
for the eigenstates using imaginary time evolution. Recently, it has been extended to 



ID lattice systems with nearest-neighbor interactions of infinite size lfT6lfT71fT8l . It has 
been successfully tested on the integrable Ising model. It is yet to be seen whether really 
difficult problems can be attacked by this algorithm. The present work provides a partial 
answer to this question by addressing the challenging problem of the quantum phase 
transition of the BHM. 

THE BOSE-HUBBARD PHASE TRANSITION 

The Hamiltonian of the BHM, in dimensionless variables, reads 

M tj M M 

H=-J^ (K+A + h - C + i E M*r - 1) - fl £ hr, (1) 

r=l 1 r=\ r=l 

where M is the number of sites. We put it finite just for the introduction, as we shall 
consider M = °° for the actual calculation, a]. (a r ) is the creation (annihilation) operator 
at site r, while h r = a\.a r . The total number of atoms N = Y!f=\ commutes with the 
Hamiltonian. In the thermodynamic limit of infinite systems, it is convenient to consider 
the particle density p = (N/M) and the chemical potential /i. 

Two parameters describe the model: the tunnelling rate J and the interaction strength 
U . Since physics of the BHM depends on the J/U ratio only, we set U = 1 for conve- 
nience. The behaviour of the system as J/U is varied has been extensively discussed 

in a number of papers, see(ISl[I2lII3EaiZDI22ll23ISl^llSll^l^llSD, thus we 
mention the most relevant details only. The Mott insulator phase occurs at integer oc- 
cupation of sites for sufficiently small tunnelling. Within the Mott insulator regime, the 
density stays constant - the system is not compressible. Typically, at the phase boundary, 
the density changes and the superfluid regime is reached via a standard quantum phase 
transition. The interesting situation occurs when the phase transition is reached via a 
path of a constant integer density p. The model belongs then to the universality class of 
the XY model ifTOl and the phase transition is of the Kosterlitz-Thouless (KT) PT1I321 
type. 

In the present paper, we determine very accurately the / value at which the KT 
transition occurs for p = 1. This is a difficult task previously addressed using dif- 
ferent methods yielding quite different results. The mean field theory [fTTl l33l yields 
J c ~ 0.086. The Bethe ansatz allows for an analytical albeit approximate treatment lfT2ll 
giving J c =l/ (2\/3) ~ 0.289. Early numerical values obtained via QMC [ 19] as well as 
renormalization group arguments [20] yielded a quite different value J c = 0.215. Later 
numerical approaches came closer to the Bethe ansatz result. In particular, Kashurnikov 
et al. flUl found J c = 0.300 ± 0.005 in a QMC study while an exact diagonaliza- 
tion for M = 12 combined with renormalization group by the same group j|22| yields 
J c = 0.304 ±0.002. Other QMC approach flU gave J c = 0.305 ±0.004. Earlier di- 
agonalization results [23] gave J c ~ 0.275 ±0.005. High order perturbative expansion 
suggests J c = 0.26 ±0.01 041 in an excellent agreement with finite DMRG treatment 
ll25Tl yielding essentially the same result. Yet other DMRG approaches gave interestingly 
different results. The periodic boundary conditions DMRG 11261 gave J c ~ 0.298 while 
the infinite size DMRG with periodic boundaries gave J c = 0.277 ± 0.01 11271 . The same 



authors attempted a finite size DMRG with open boundary conditions (yielding a better 
convergence) giving J c = 0.297 ±0.01 ll28l . Therefore, the current algorithms disagree 
on the the value of J c . This is partially due to logarithmic finite-size effects close to the 
critical point of the KT transition 1)221 l28l . The accurate determination of J c is thus an 
ideal testbed for algorithms dealing with infinite systems. 



TIME-EVOLVING BLOCK DECIMATION ALGORITHM 

The approach |[T6ll used is an infinite variant of the TEBD algorithm lfT4l[l"5ll . It belongs 
to a rapidly developing class of methods which utilise the understanding of quantum en- 
tanglement in complex systems to facilitate effective simulations of many-body quantum 
systems. All these methods have some links (see discussion in 11341 ) with the DMRG al- 
gorithms [13] extending them either to time evolution [TT4l[T5l and to many dimensional 
systems 0511361 1371 ) A particularly interesting variant of the TEBD approach assumes 
translational invariance of the wavefunction and enables a treatment of explicitly infinite 
systems lfT6l[T71[T8ll . While we refer the reader to original papers for details, we sketch 
here the main idea. 

The BHM Hamiltonian involves only nearest-neighbor interactions and can thus be 
written as a sum over consecutive sites r of reduced hamiltonians M r,r+1 ). Each site is 
described by a Hilbert space of dimension D. We follow the imaginary time evolution, 
i.e., compute (we assume fi = 1 in the whole paper) 

= ap(-Ht)|<P(0)> 
1 " ||exp(-Hl)|W(0)>|| 1 ; 

For large % and a generic initial guess ^(O)), such a procedure should yield a good 
approximation of the ground state provided there is a gap in the spectrum of the system. 

We assume translational invariance of our system, i.e., /z( r ' r+1 )'s does not depend on 
site r. Moreover we restrict ourselves to wavefunctions that are invariant under shifts by 
one lattice spacing^] The time evolution proceeds in the following way. The normalized 
wave function is Schmidt decomposed at site r as 

m=t^ r) \* { a r) )®\*a r+ % (3) 
ce=l 

where % is a finite Schmidt rank, the |<J> a ) states are an orthonormal set and Ea(^«) 2 = 1 
for global normalization. As discussed in [[141 1X311 . for low lying states, the X a values 
decrease rapidly, so that a good approximation of the state is obtained keeping only the 
few largest A„. The two consecutive left or right sub lattices states are related via tensors 



1 This assumption is not essential for the method and can be easily lifted, as pointed out in [ 16]. It excludes 
"charge density wave" or "checkerboard" solutions. Previous studies of the BHM indicate that the ground 
state is translationally invariant so there is no real limitation in the approach used. 



r( r ) as follows: 



j8=n=i 

j3=li=l 

where \i^ r ') is an orthonormal basis spanning the local Hilbert space at site r (e.g. the 
Fock basis for the BHM). In that way one obtains a so-called matrix product state (MPS) 
ll38Tl representation of |*P(t)). The imaginary time evolution operation can be expressed 
as a product of elementary two-site evolution operators 

V {r,r+\) =exp (-h^ r + l )8t), (5) 

using the Suzuki-Trotter decomposition [|39l . 

Now the crucial observation is that, for translationally invariant \*¥), all and X^ 
are independent of r. The application of a two-site gate [/( r,r+1 ) has two effects: 

• it breaks for a moment the translational invariance (since sites r and r + 1 are 
updated only); 

• the resulting vector is not longer a MPS. 

The latter is taken care of by a singular value decomposition of the resulting matrix 
and retaining only the % largest eigenvalues (and corresponding vectors) - see |[T4l [T5l 
for details and lfT6l for a graphical representation. This approximation brings the vector 
back into the MPS form. The modified T^'s, A^'s are then updated for all the sites - 
the translational invariance is restored. 

It is worth stressing that the repetition of this procedure allows one to act on two 
consecutive sites only, taking proper care of local updates of appropriate tensors. In 
the limit % — > °° and 8t — > 0, the exact evolution is obtained. The error introduced 
by finiteness of 8t is easily controlled (it varies as a power of 8t depending on the 
order of the Suzuki-Trotter decomposition) while the strength of the method is that the 
number % of not too small X values remains reasonably small. As shown in IfToll . to 
facilitate convergence, one can take successively smaller and smaller time steps. In the 
following we use a simple second order Trotter expansion (which, according to our 
experience, seems to be the most efficient in imaginary time evolution when time steps 
are decreased). 

APLICATION TO THE BOSE-HUBBARD MODEL 

To make the infinite-TEBD method work, additional parameters have to be adjusted. One 
of them is the maximal allowed occupation N max at a given site (giving the dimension 
of the Hilbert space at each site, D = N max + 1). In earlier DMRG calculations for mean 
unit occupancy typically N max = 4 or N max = 6 was used [|25l |28]|. We take N max = 6, 
although the difference with N max = 4 is small. The second important cutoff parameter 



is the restriction on %, the number of A values retained at each step. This value can be 
quite small deep in the Mott insulator region, it has to be taken painfully large close to the 
KT transitions. This is costly since the execution time scales as % 3 IfToll . A satisfactory 
convergence of the correlation function 

C(s) =< al +r a r > (6) 

requires % of the order of one hundred. Such correlation functions define the phase 
coherence of the condensate, we shall use them to determine the transition point. 

In the Mott regime, the system has a finite gap, which implies that |^P(t)), eq. ([5]), con- 
verges exponentially fast at long times to the ground state, at a rate proportional to the 
energy gap. This is indeed what we numerically observe, with exponential convergence 
of the energy too. Thus, gor a given time step we perform the imaginary time evolution 
until the slope of the energy curve w.r.t. imaginary time reaches a desired small thresh- 
old. Then the time step is decreased by a factor 2 and the procedure repeated until a 
computer accuracy limit for the time step is reached. While this procedure may not be 
optimal for speed, it assures the convergence of the results obtained. 

It seems quite surprising that our results both for the energy and for the correlation 
function, Fig. [TJ converge well even in the superfluid regime where the gap in the 
excitation spectrum should vanish. We indeed observe such a convergence, although 
apparently not fully exponential. We presently have no clear explanation of this behavior, 
which is under study. A plausible cause is the enforcement of translational invariance at 
each step of the algorithm. That restricts the spanned Hilbert space to the translationally 
invariant subspace killing, e.g., the possible excitations of low-frequency Bogolyubov 
modes in the superfluid regime. We carefully checked that all results presented in this 
paper are converged with respect to a modification of the (small) time step and an 
increase of %. 

The convergence of the long-range correlation functions is of particular importance. 
As discussed in E71l28ll . at integer density, it should exhibit a power law C(s) s~ K ^ 2 for 
large s, on the superfluid side with K increasing up to K = 0.5 when the KT transition is 
approached. This is a behaviour typical for a Luttinger liquid characterising low energy 
excitations in the superfluid state. Fig. [2] reveals that it is indeed the case. As expected, 
the power law behaviour no longer holds on the Mott side of the transition. Here, deep 
in the Mott regime, an exponential decay of correlations is observed. 

To estimate the KT transition point, a special case must be paid to keep a constant 
unit density. While for a finite system the particle number conservation can be explicitly 
used to keep the density fixed this is not the case for the infinite version of the code IfToll . 
Instead, a chemical potential, /i, must be used and adjusted to keep a contant desided 
density. Our simulations were performed for different values of \x and the correlation 
functions shown (and used later to extract the transition point) were obtained for p well 
within [0.999, 1.001] interval. To obtain slopes at p = 1 an appropriate interpolation of 
the results has been used. We observed that the power law coefficient in that interval 
changed weakly, the possible error in the transition point evaluation is included in the 
estimates given below. 

The correlation functions C[s) obtained at different J values are fitted with a power 
law on the superfluid side of the transition point. Since in the Mott regime the behaviour 




FIGURE 1. Correlation function C(s) in the superfiuid phase regime, J = 0.32. Solid line % = 30, 
dashed % = 100 while the dots and the underlying thin solid line correspond to % — 150 and % = 200, 
respectively. The last two sets almost coincide indicating the convergence with respect to %. The linear 
behaviour in this log-log plot proves that the correlation function decays algebraically. 



of the correlation function hardly resembles a power law (see Fig. [2]), we make fits in 
different ranges of s exactly following the procedure of Kiiehner et al [28 J. The obtained 
slopes (as a function of J) are interpolated to yield the transition point J c . These points 
are collected in Table I for different s intervals. For comparison we reproduce also the 
results of [|28]. While earlier DMRG results seem to be quite sensitive to the range of 
s values, our fits yield much closer and mutually consistent results. The comparison 
presented shows the power and accuracy of the infinite TEBD algorithm. 

Clearly, fits for small s lead to misleading results, C(s) for s below few tens (say 30) 
does not yet behave in the asymptotic power law way even in the superfiuid regime. 
Given the accuracy of the method we can extend the evaluation of C(s) to much larger 
s values, yielding global fits reported at the bottom of Table I. Note that those fits are in 
decent agreement with moderate s values of the order of 50. 

We cannnot, however, extend the calculation of C(s) to much larger s values of the 
order of several thousands or more due to the numerical precision. 

The KT transition point at unit density should be safely close to the value given by 
the last fit J c = 0.2975 ± 0.0005 where the (conservative) estimate of the error is coming 
from the comparison with other fits. While the obtained value is quite close to that 
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FIGURE 2. Correlation functions C(s) for different values of J on a double logarithmic scale. From 
bottom to top: J = 0.26,0.29,0.293,0.296,0.30. The J = 0.26 curve in the deep Mott regime regime 
corresponds to a Mott insulator. As expected it resembles an exponential. In the deep superfluid regime, 
the power law is very well obeyed, that allows for accurate fits. Close to the transition, even on the Mott 
side, C(s) shows transient behaviour for small s, then a region of an approximate power law behaviour, 
then for very large s an exponential tail is expected in the Mott regime. When approaching the critical 
point on the Mott side, this exponential region takes place at larger and larger s. In the range s < 800 
shown in this plot, it is not visible expect for the lowest / = 0.26 curve. The J = 0.296 curve, the closest 
one to the critical curve s~'/ 4 , is shown by a dashed line. 

reported previously (281, even our conservative error estimate makes the result more 
accurate by one order of magnitude. 

Quite a similar approach can be used for a double occupancy case p = 2. Here, using 
the very same method, we obtain a benchmark value / c (2) = 0.175 ±0.002. For that 
calculation the possible maximal occupancy of each site had to be larger or equal to 6. 
We have checked by a comparison of selected runs with those for even larger occupancy 
that the assumed maximal occupation is sufficient for the convergence. 

CONCLUSIONS 

To summarize, we have shown that the infinite TEBD algorithm lfl6l may be efficiently 
applied to nontrivial challenging problems such as the determination of the KT critical 



TABLE 1. The location of the critical point J c , esti- 
mated from different ranges of s value in the correlation 
function C(s). The second column reproduces fits of [28], 
the third column yields present results. 

s J c (DMRG) J c 

4<s<8 0.2874 ±0.0001 0.2868 ±0.0006 

8 < s < 16 0.2938 ± 0.0001 0.2933 ± 0.0002 

16 < .v < 32 0.2968 ±0.0003 0.2955 ±0.0001 

32<5<48 0.3062±0.0003 0.2970±0.0001 

48<s<64 0.3 107 ±0.01 0.2975 ±0.0001 

100<i<250 0.2980±0.0001 

100 < s < 500 0.2972 ± 0.0001 

50<s<500 0.2975 ±0.0001 



point in the BHM. The accuracy is superior to that obtained by the standard DMRG 
approach. Surprisingly, converged results are also obtained in the superfluid regime. 
The restriction of the algorithm to the translationally invariant subspace, while justified 
for present ground state studies, may restrict the application of the method ifTBI to the 
real time evolution. There it could allow to treat "exactly" the dynamics of the quantum 
phase transition. Work in this direction is in progress. 
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